# scp savepoints/savepoint_2023-05-15/controls.rds UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/.
# scp savepoints/savepoint_2023-05-15/ww1.rds UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/.
# scp UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/ma5.3.1.rds savepoints/savepoint_2023-05-15/.
# scp UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/ma5.3.2.rds savepoints/savepoint_2023-05-15/.
# scp UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/ma5.3.3.rds savepoints/savepoint_2023-05-15/.
# scp UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/ma5.4.1.rds savepoints/savepoint_2023-05-15/.
if(!exists("controls")) controls = readRDS(fs::path("../savepoints/savepoint_2023-05-15/controls.rds"))
source("setup.R")
ww1 = readRDS(fs::path("../",controls$savepoint,"ww1.rds"))
shapes = readRDS(fs::path("../",controls$savepoint,"shapes.rds"))We now model all ARAs together.
ww_all = ww1 %>%
# log
dplyr::mutate(logvl=log(vl)) %>%
mutate(vl=if_else(vl==0, 1, vl)) %>%
# create indexes for INLA
dplyr::mutate(day1=day,
ara1=as.numeric(as.factor(ara_n)),
ara2=ara1) %>%
# group KLBS and KLZH (otherwise not identifiable)
dplyr::mutate(lab2=if_else(lab=="KLBS","KLZH",lab),
lab_n2=as.numeric(as.factor(lab2))) %>%
# group lab and method
dplyr::mutate(lab_method=paste0(lab2,"_",method),
lab_method_n=as.numeric(as.factor(lab_method)))
# correspondence table
corr_all = ww_all %>%
group_by(ara_n,ara_id,ara_name,kt,pop,lab,lab_n,lab2,lab_n2,lab_method,lab_method_n,ara1,ara2,NUTS2_name) %>%
count() %>%
ungroup()
corr_all_ara = ww_all %>%
group_by(ara1,ara_name,ara_id,kt,NUTS2_name) %>%
count() %>%
ungroup()
if(!controls$rerun_models) {
mw_100_desc_table(ww_all) %>%
dplyr::mutate(across(everything(),as.character)) %>%
tidyr::gather() %>%
flextable::flextable(cwidth=c(4,4))
}key | value |
|---|---|
Number of ARAs | 118 |
Number of laboratories | 9 |
Number of laboratory methods | 11 |
Number of measurements | 20535 |
Measurements below LOQ | 235 |
Measurements below LOD | 97 |
First | 2022-02-07 |
Last | 2023-05-14 |
Median viral concentration [gc/L] | 1e+05 (range: 0 to 8e+06) |
Median flow [m3/day] | 1e+04 (range: 3e+02 to 9e+05) |
Median viral load [gc/day/100,000] | 4e+12 (range: 1 to 7e+14) |
We directly apply model A5.3 to all ARAs. We assume one temporal trend for Switzerland, and each ARA is free to deviate from it independently.
if(controls$rerun_models) {
ma5.3.1 = INLA::inla(vl ~ 1 +
f(below_loq,model="iid") +
f(below_lod,model="iid") +
f(day,model="rw2", scale.model=TRUE, constr=TRUE,
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
f(hol,model="linear",mean.linear=0,prec.linear=.2) +
f(ara1,model="iid") +
f(day1,model="rw1", scale.model=TRUE, constr=TRUE,
group=ara2, control.group=list(model="iid"),
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))),
data = ww_all,
family = "gamma",
control.compute = list(waic=TRUE,config=TRUE),
control.predictor = list(compute=TRUE,link=1))
saveRDS(ma5.3.1,file=paste0("../",controls$savepoint,"ma5.3.1.rds"))
} else {
ma5.3.1 = readRDS(file=paste0("../",controls$savepoint,"ma5.3.1.rds"))
}
summary(ma5.3.1)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", "
## blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)")
## Time used:
## Pre = 1.08, Running = 12194, Post = 12.8, Total = 12208
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 13.925 7.710 -1.447 13.925 29.296 13.926 0
## weekend -0.002 0.009 -0.020 -0.002 0.016 -0.002 0
## hol -0.001 0.026 -0.053 -0.001 0.050 -0.001 0
##
## Random effects:
## Name Model
## below_loq IID model
## below_lod IID model
## day RW2 model
## ara1 IID model
## day1 RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 3.448 0.048 3.351 3.450 3.538 3.456
## Precision for below_loq 3.111 1.421 1.129 2.852 6.613 2.367
## Precision for below_lod 0.012 0.006 0.004 0.011 0.029 0.009
## Precision for day 0.011 0.003 0.008 0.011 0.018 0.010
## Precision for ara1 1.064 0.253 0.702 1.017 1.683 0.916
## Precision for day1 0.482 0.020 0.442 0.483 0.521 0.485
##
## Watanabe-Akaike information criterion (WAIC) ...: 1223613.19
## Effective number of parameters .................: 4231.92
##
## Marginal log-Likelihood: -732342.88
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## exp(beta) 0.025quant 0.975quant
## weekend 1 0.98 1.02
## hol 1 0.95 1.05
We add a covariate to measure the effect of the lab and of changes in methodology. To allow identifiability, we have to make sure that there are multiple ARAs per laboratory, so we group together KLBS (only 1 ARA) and KLZH (12 ARAs). The reference lab is ALTGR.
if(controls$rerun_models) {
ma5.3.2 = INLA::inla(vl ~ 1 +
f(below_loq,model="iid") +
f(below_lod,model="iid") +
f(day,model="rw2", scale.model=TRUE, constr=TRUE,
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
f(hol,model="linear",mean.linear=0,prec.linear=.2) +
f(ara1,model="iid") +
f(day1,model="rw1", scale.model=TRUE, constr=TRUE,
group=ara2, control.group=list(model="iid"),
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
lab_method,
data = ww_all,
family = "gamma",
control.compute = list(waic=TRUE,config=TRUE),
control.predictor = list(compute=TRUE,link=1))
saveRDS(ma5.3.2,file=paste0("../",controls$savepoint,"ma5.3.2.rds"))
} else {
ma5.3.2 = readRDS(file=paste0("../",controls$savepoint,"ma5.3.2.rds"))
}
summary(ma5.3.2)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", "
## blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)")
## Time used:
## Pre = 4.71, Running = 16149, Post = 4.68, Total = 16159
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 14.624 3015.142 -5898.535 14.624 5927.783 14.624 0
## lab_methodEAWAG_0 0.467 0.755 -1.012 0.467 1.947 0.467 0
## lab_methodEurofins_0 -3.076 0.574 -4.201 -3.076 -1.951 -3.076 0
## lab_methodEurofins_1 -1.522 0.577 -2.653 -1.522 -0.390 -1.522 0
## lab_methodKLZH_0 0.698 0.615 -0.509 0.698 1.904 0.698 0
## lab_methodLdU_0 0.019 0.654 -1.264 0.019 1.302 0.019 0
## lab_methodMicrosynth_0 -0.892 0.504 -1.880 -0.892 0.096 -0.892 0
## lab_methodMicrosynth_1 -1.262 1.694 -4.585 -1.262 2.061 -1.262 0
## lab_methodSCAV_NE_0 -1.125 0.980 -3.046 -1.125 0.797 -1.125 0
## lab_methodSUPSI_0 -0.713 0.809 -2.300 -0.713 0.874 -0.713 0
## weekend -0.004 0.009 -0.022 -0.004 0.013 -0.004 0
## hol -0.007 0.025 -0.057 -0.007 0.042 -0.007 0
##
## Random effects:
## Name Model
## below_loq IID model
## below_lod IID model
## day RW2 model
## ara1 IID model
## day1 RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 3.290 0.033 3.232 3.288 3.362 3.279
## Precision for below_loq 1.579 0.553 0.601 1.519 2.690 1.448
## Precision for below_lod 0.000 0.001 0.000 0.000 0.000 0.000
## Precision for day 0.131 0.012 0.105 0.131 0.151 0.137
## Precision for ara1 0.263 0.099 0.154 0.240 0.528 0.187
## Precision for day1 0.512 0.017 0.477 0.513 0.544 0.515
##
## Watanabe-Akaike information criterion (WAIC) ...: 1223815.85
## Effective number of parameters .................: 4045.45
##
## Marginal log-Likelihood: -732451.58
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## exp(beta) 0.025quant 0.975quant
## lab_methodEAWAG_0 1.60 0.36 7.01
## lab_methodEurofins_0 0.05 0.01 0.14
## lab_methodEurofins_1 0.22 0.07 0.68
## lab_methodKLZH_0 2.01 0.60 6.71
## lab_methodLdU_0 1.02 0.28 3.67
## lab_methodMicrosynth_0 0.41 0.15 1.10
## lab_methodMicrosynth_1 0.28 0.01 7.85
## lab_methodSCAV_NE_0 0.32 0.05 2.22
## lab_methodSUPSI_0 0.49 0.10 2.40
## weekend 1.00 0.98 1.01
## hol 0.99 0.94 1.04
Instead of just one temporal trend for Switzerland, we now allow independent temporal trends for each NUTS-2 region. ARAs- within each region are then allowed to deviate from the regional trend independently.
if(controls$rerun_models) {
ma5.3.3 = INLA::inla(vl ~ 1 +
f(below_loq,model="iid") +
f(below_lod,model="iid") +
f(day,model="rw2", scale.model=TRUE, constr=TRUE,
group=NUTS2, control.group=list(model="iid"),
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
f(hol,model="linear",mean.linear=0,prec.linear=.2) +
f(ara1,model="iid") +
f(day1,model="rw1", scale.model=TRUE, constr=TRUE,
group=ara2, control.group=list(model="iid"),
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
lab_method,
data = ww_all,
family = "gamma",
control.compute = list(waic=TRUE,config=TRUE),
control.predictor = list(compute=TRUE,link=1))
saveRDS(ma5.3.3,file=paste0("../",controls$savepoint,"ma5.3.3.rds"))
} else {
ma5.3.3 = readRDS(file=paste0("../",controls$savepoint,"ma5.3.3.rds"))
}
summary(ma5.3.3)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", "
## blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)")
## Time used:
## Pre = 3.02, Running = 10283, Post = 4.05, Total = 10290
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 14.766 16.104 -19.590 14.766 49.118 14.768 0
## lab_methodEAWAG_0 0.358 0.188 -0.010 0.358 0.730 0.357 0
## lab_methodEurofins_0 -3.035 0.181 -3.389 -3.035 -2.678 -3.036 0
## lab_methodEurofins_1 -1.838 0.196 -2.222 -1.839 -1.452 -1.839 0
## lab_methodKLZH_0 0.646 0.173 0.307 0.645 0.987 0.644 0
## lab_methodLdU_0 -0.148 0.175 -0.489 -0.149 0.197 -0.151 0
## lab_methodMicrosynth_0 -1.031 0.138 -1.300 -1.032 -0.760 -1.032 0
## lab_methodMicrosynth_1 -1.376 0.644 -2.639 -1.376 -0.113 -1.376 0
## lab_methodSCAV_NE_0 -1.078 0.265 -1.599 -1.078 -0.558 -1.077 0
## lab_methodSUPSI_0 -0.933 0.224 -1.372 -0.933 -0.492 -0.934 0
## weekend -0.003 0.009 -0.020 -0.003 0.015 -0.003 0
## hol 0.009 0.026 -0.042 0.009 0.059 0.009 0
##
## Random effects:
## Name Model
## below_loq IID model
## below_lod IID model
## day RW2 model
## ara1 IID model
## day1 RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 3.301 0.043 3.223 3.298 3.391 3.290
## Precision for below_loq 0.780 1.139 0.076 0.447 3.574 0.185
## Precision for below_lod 0.004 0.004 0.000 0.002 0.015 0.001
## Precision for day 0.005 0.000 0.004 0.005 0.006 0.005
## Precision for ara1 9.554 1.243 7.164 9.537 12.048 9.630
## Precision for day1 0.951 0.052 0.858 0.948 1.062 0.939
##
## Watanabe-Akaike information criterion (WAIC) ...: 1223915.54
## Effective number of parameters .................: 3677.11
##
## Marginal log-Likelihood: -748599.68
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## exp(beta) 0.025quant 0.975quant
## lab_methodEAWAG_0 1.43 0.99 2.07
## lab_methodEurofins_0 0.05 0.03 0.07
## lab_methodEurofins_1 0.16 0.11 0.23
## lab_methodKLZH_0 1.91 1.36 2.68
## lab_methodLdU_0 0.86 0.61 1.22
## lab_methodMicrosynth_0 0.36 0.27 0.47
## lab_methodMicrosynth_1 0.25 0.07 0.89
## lab_methodSCAV_NE_0 0.34 0.20 0.57
## lab_methodSUPSI_0 0.39 0.25 0.61
## weekend 1.00 0.98 1.01
## hol 1.01 0.96 1.06
#’ ## Model A5.4.1: spatial structure #’ #’ #+ ma5.4.1a, fig.width=8, fig.height=12, R.options = list(width = 1000) if(controls\(rerun_models) { shapes_all = shapes\)ara_shp %>% dplyr::filter(ara_id %in% ww_all\(ara_id) %>% left_join(corr_all_ara,by = join_by(ara_id)) %>% dplyr::arrange(ara1) sf_use_s2(FALSE) graph_all = spdep::poly2nb(shapes_all) path_graph = paste0("../",controls\)savepoint,“W_all.adj”) nb2INLA(path_graph, graph_all) ma5.4.1 = INLA::inla(vl ~ 1 + f(below_loq,model=“iid”) + f(below_lod,model=“iid”) + f(day,model=“rw2”, scale.model=TRUE, constr=TRUE, group=NUTS2, control.group=list(model=“iid”), hyper=list(prec = list(prior = “pc.prec”, param = c(1, 0.01)))) + f(weekend,model=“linear”,mean.linear=0,prec.linear=.2) + f(hol,model=“linear”,mean.linear=0,prec.linear=.2) + f(ara1,model=“bym2”, graph=path_graph, scale.model = TRUE, constr = TRUE, hyper = list(theta1 = list(“PCprior”, c(1, 0.01)), theta2 = list(“PCprior”, c(0.5, 0.5)))) + f(day1,model=“rw1”, scale.model=TRUE, constr=TRUE, group=ara2, control.group=list(model=“iid”), hyper=list(prec = list(prior = “pc.prec”, param = c(1, 0.01)))) + lab_method, data = ww_all, family = “gamma”, control.compute = list(waic=TRUE,config=TRUE), control.predictor = list(compute=TRUE,link=1)) saveRDS(ma5.4.1,file=paste0(“../”,controls\(savepoint,"ma5.4.1.rds")) } else { ma5.4.1 = readRDS(file=paste0("../",controls\)savepoint,“ma5.4.1.rds”)) } summary(ma5.4.1) summary_exp_vl(ma5.4.1,pars=“lab|method|hol|weekend”) ppp_vl_ara(ww_all,ma5.4.1) #+ ma5.4.1b, fig.width=8, fig.height=6, R.options = list(width = 1000) mw_130_map_relative_vl(ma5.4.1,corr_all_ara,shapes) #+ ma5.4.1c, fig.width=8, fig.height=6, R.options = list(width = 1000) mw_131_map_deviation_from_average(ma5.4.1,corr_all_ara,ww_all,shapes,10)